#ifndef BIO_PH_CHEMSYS_INCLUDED
#define BIO_PH_CHEMSYS_INCLUDED 1

#include "kernel/reaction.hpp"
#include "yocto/math/kernel/matrix.hpp"
#include "yocto/math/fcn/derivative.hpp"
#include "yocto/math/kernel/linsys.hpp"

typedef set<string,Reaction> ReactionSet;

class Initializer;
class Solution;

class ChemSys
{
public:
	explicit ChemSys( const Library &lib );
	virtual ~ChemSys() throw();
	
	void append( const Reaction &r );
	void load( lua_State *L, const string &chemsys_name );
	
	void compile();
	const size_t N; //!< #reactions
	const size_t M; //!< library.size() 
	
	void randomX() throw();
	void computeGammaAndPhi(Real t, bool compute_dtGam); //!< From Current X => Gamma, Phi
	void computeXi();                //!< from Phi and Nu => Nu'* (Phi*Nu')^{-1}
	void computeChi();               //!< from Xi         => Nu'* (Phi*Nu')^{-1} * Phi
	
	void solveGamma( Real t );       //!< Newton step
	
	void updateRate( Solution &dS, Real t , const Solution &S);
	
	void load_data( const array<Real> &F ); //!< F.size() >= M
	void save_data( array<Real>  &F) const; //!< F.size() >= M;
	
	const Reaction & operator[]( const string &id ) const;
	const Reaction & operator[]( const char   *id ) const;
	
	
private:
	YOCTO_DISABLE_COPY_AND_ASSIGN(ChemSys);
	ReactionSet   reactions_;
	
public:
	const Library &library;
	
private:
	vector<Function> Kproc; //!< function to build K, size=N
	vector<Real>     X;     //!< vector of concentrations, size=M
	vector<Real>     Gamma; //!< the Gamma value (size=N)
	vector<Real>     dtGam; //!< dGamma/dt (size=N)
	matrix<Real>     Nu;    //!< algebraic stochiometic coefficients (size=NxM)
	matrix<Real>     NuP;   //!< product stochiometric  (size=NxM)
	matrix<Real>     NuR;   //!< -product stochiometrix (size=NxM)
	matrix<Real>     Phi;   //!< dGamma/dX (size=NxM)
	matrix<Real>     W;     //!< Phi * Nu' (size=NxN), then LU
	matrix<Real>     iW;    //!< inverse of (Phi*Nu')
	linsys<Real>     lss;   //!< to perform decomposition
	matrix<Real>     Xi;    //!< Nu'* (Phi*Nu')^{-1} (size=MxN)
	matrix<Real>     Chi;   //!< Nu'* (Phi*Nu')^{-1} * Phi (size=MxM)
	vector<Real>     dX;    //!< for Newton step
	
	friend class Initializer;
	friend class Solution;
	
public:
	derivative<Real> drvs;  //!< derivator
	Real             dt;    //!< for derivative of constants w.r.t. time.
	Real             ftol;  //!< fractional tolerance for a full Newton step and other algorithm
	
	inline const array<Real>  & _X()     const throw() { return X; }
	inline const array<Real>  & _Gamma() const throw() { return Gamma; }
	inline const matrix<Real> & _Xi()    const throw() { return Xi; }
	inline const matrix<Real> & _Phi()   const throw() { return Phi; }
	inline const array<Real>  & _dtGam() const throw() { return dtGam; }
};


#endif
